Culture and multiomic analysis of lung cancer patient-derived pleural effusions revealed distinct druggable molecular types

Malignant pleural effusion (MPE) is an independent determinant of poor prognostic factor of non-small cell lung cancer (NSCLC). The course of anchorage independent growth within the pleural cavity likely reforms the innate molecular characteristics of malignant cells, which largely accounts for resistance to chemotherapy and poor prognosis after the surgical resection. Nevertheless, the genetic and transcriptomic features with respect to various drug responses of MPE-complicated NSCLC remain poorly understood. To obtain a clearer overview of the MPE-complicated NSCLC, we established 28 MPE-derived lung cancer cell lines which were subjected to genomic, transcriptomic and pharmacological analysis. Our results demonstrated MPE-derived NSCLC cell lines recapitulated representative driver mutations generally found in the primary NSCLC. It also exhibited the presence of distinct translational subtypes in accordance with the mutational profiles. The drug responses of several targeted chemotherapies accords with both genomic and transcriptomic characteristics of MPE-derived NSCLC cell lines. Our data also suggest that the impending drawback of mutation-based clinical diagnosis in evaluating MPE-complicated NSCLS patient responses. As a potential solution, our work showed the importance of comprehending transcriptomic characteristics in order to defy potential drug resistance caused by MPE.

Lung cancer is one of the most prevalent cause of cancer-related death worldwide, with an estimated one million deaths annually. Within Korean population, it is the most lethal disease and the third common types of cancer 1,2 . It is also the most prevailing cause of malignant pleural effusion (MPE), accounting for approximately 30% of such effusions 3 . Although recent advances in prevention and early detection have contributed to the general increase in the survival rate, the occurrence of MPE in the context of lung cancer designates terminal stage and an independent prognostic parameter of the poorer survival rate 4 . Lung cancer patients with MPE mostly exhibit insensitivities to certain chemotherapies mainly due to high tumor burden 5,6 . Besides, the course of acquiring anoikis resistance alters multiple tumorigenic pathways, which eventually brings about insensitivity to certain chemotherapies 7 . Nevertheless, the genetic and transcriptomic features with respect to various drug responses of MPE-complicated NSCLC have not been thoroughly studied yet, which makes the implication of drug insensitivities of MPE-derived tumor cell still hypothetical.
Moreover, several studies have indicated that acquired resistance to second line chemotherapy such as crizotinib is associated with aggressive disease progressions 8 . Nevertheless, the widely accepted oncogenic association of secondary or tertiary alternations by the second line chemotherapy has been focused on the primary NSCLC, and oncogenic impact of chemotherapy-derived molecular changes in MPE remains unclear.
To obtain a clearer overview of the MPE-complicated NSCLC in accordance with heterogeneous drug responses, we established 28 lung cancer cell lines derived directly from the MPE, and two crizotinib-resistant sublines which were then subjected to genomic, transcriptomic and pharmacological analysis. Our results not only showed that MPE-derived NSCLC cell lines recapitulated representative genetic aberrations of the primary MPE-derived lung cancer cell lines retain the representative genetic mutations of human lung adenocarcinoma. Several studies have shown that MPE-derived tumor cell lines maintain the mutational profiles of the lung adenocarcinoma 6,9 . We performed both whole exome sequencing (WES) and targeted sequencing (CCP) on the 28 MPE-derived cell lines. The overall mutational spectrum of the MPE-derived NSCLC cell lines such as mutant-allele tumor heterogeneity (MATH) analysis and relative contribution of mutational signatures was analyzed using WES data. The targeted sequencing as well as direct Sanger sequencing re-confirmed the presence of key driver mutations minimizing the possible false-positive sequencing errors ( Supplementary Fig. 2, Supplementary Table 2). The precise exclusion of the germinant mutations was not feasible since the DNA from matched blood or normal tissue was unavailable when WES of established cell lines was performed. Alternatively, we have referred to the Clinvar database (https:// www. ncbi. nlm. nih. gov/ clinv ar) in order to avoid misestimating the mutation frequency due to inflated germinant/benign mutations. The pathogenicity of the representative driver mutations was manually inspected, and all marked mutations for calculating the frequency of mutated tumor genes in Fig. 2A were previously reported as "Pathogenic", "Likely Pathogenic" or "Drug Response" in the Clinvar database.
To establish a more comprehensive mutational profile, the aberrations in key driver genes were compared to a larger lung cancer cell line cohort from the cancer cell line encyclopedia project (CCLE) 10 which was composed of a total of 63 lung cancer cell lines ( Fig. 2A). As a first validation, our cohort had a significantly higher rate of EGFR mutations (36%) compared to CCLE cohort (14%). It has been reported that the EGFR mutation frequency in lung adenocarcinoma is approximately 10% as much larger cohort (TCGA) 11 indicated ( Supplementary Fig. 3). Pervasive EGFR mutation in our cohort was accordant with the previous MPE of lung adenocarcinoma study 12 . In contrast, KRAS mutation was detected at a low rate in our cohort (7%) compared to CCLE (35%) and TCGA (30%) cohorts, which was unexpected since mutant KRAS is reported to promote malignant pleural effusion formation 5 . Previous study indicated that EGFR is more frequently aberrated than KRAS in Korean NSCLC population 13 , which suggest that the geographic factor partially shaped the mutational landscape of our cohort. Nine patients harbored driver fusion genes including EML4-ALK, KIF5B-RET and CD74-ROS1 in our cohort ( Supplementary Fig. 4). Few groups reported that lung cancer patients with fusion genes, especially EML4-ALK were younger than those without 14,15 , which was in parallel with our data. The presence of fusion genes was mutually exclusive with EGFR and KRAS mutation in parallel with prior study 14 .
We also calculated the heterogeneity scores by using a previously reported method, mutant-allele tumor heterogeneity (MATH) 16 . This analysis allowed to estimate the dispersion of mutant allele frequencies (MAFs) within a cell population, which reflected the intra-tumor genetic heterogeneity. The cell lines with mutant EGFR had higher MATH score than two other groups, suggesting that higher level of genetic heterogeneity is causative factor for EGFR mutation (Fig. 2B). Then, we evaluated the relative contribution of single nucleotide variations (SNVs) to the mutational signature of each cell line using web-based algorithm, MUSICA. The most predominant point mutation type in the cell lines was the C-to-T transitions (Fig. 2C, Supplementary Table 3). Except for SNU-3023 that exhibited distinct configurations of the relative contribution of SNVs, all cell lines had analogous patterns. We then compared the most various mutational contribution among the three groups. A[T > C]C was the most different SNV signature followed by G[C > T]T (Fig. 2D).
Overall, comprehensive mutational analysis indicated that there is a distinct separation of genomic drivers in our cohort: EGFR aberration, fusion genes, and others. We have applied these groups for further analysis of transcriptome and drug screening. www.nature.com/scientificreports/ analysis also displayed that cell lines with the fusion genes were clustered distinctly from the cell lines with EGFR mutation (Fig. 3B). The other group consisting of cell lines with genomic aberrations other than EGFR and fusion genes were not clustered at all. Next, we performed gene set enrichment assay (GSEA) to identify pathways that are specifically enriched in each group. Only spermatogenesis pathway was significantly enriched in EGFR group compared to other two groups (p < 0.05). Multiple tumor driver genes that are associated with spermatogenesis pathway such as mTOR, EZH2, NF2, DCC and MLF1 had high enrichment score in the EGFR group (Fig. 3C). We further validated the enriched-pathways in EGFR-mutant cell lines using CCLE dataset. We performed gene set enrichment analysis (GSEA) using the CCLE dataset restricted to EGFR mutated NSCLC cell lines that were introduced in mutational analysis section ( Fig. 2A). The EGFR-mutant cell lines from the CCLE dataset were HCC2279, HCC4006, HCC827, NCIH1355, NCIH1793, NCIH1975, NCIH2291, NCIH3255 and PC14. We unified all running parameters of GSEA and the input matrix by accounting only shared genes between CCLE and our mRNA expression dataset which consists of 15,970 genes. Then, the gene sets that are upregulated in EGFR-mutated cell lines from CCLE were estimated in comparison with cell lines without EGFR mutation from our cohort. Indeed, spermatogenesis pathway was highly enriched at nominal p-value < 0.005. In contrast to our cohort, the GSEA using CCLE dataset identified G2M_checkpoint and E2F_targets pathways as significantly enriched at nominal p-value < 0.01 as well. This analysis reiterates the EGFR-mutated NSCLC indeed has distinctive transcriptomic subtype represented by spermatogenesis pathway in GSEA. Meanwhile, several pathways were enriched in the fusion gene group (Fig. 3D): Inflammatory response (p < 0.01), Hypoxia (p < 0.01), TNFα signaling (p < 0.05), Angiogenesis (p < 0.05). Genes with high enrichment scores for each sign- Induced resistance to crizotinib causes nuclear localization of β-catenin in MPE-derived NSCLC through increased DKK1 expression. An ALK inhibitor is the standard treatment for advanced NSCLC patients harboring the anaplastic lymphoma kinase (ALK) fusion gene. However, secondary ALK mutations or alternative pathway changes cause a fraction of the tumors to be resistant to ALK inhibitors 17,18 . We further validated that the acquired resistance to crizotinib, the first generation of the ALK inhibitor, alters the transcriptional patterns of MPE-derived NSCLC cell lines harboring EML4-ALK fusion gene. We established two crizotinib-resistant sublines (SNU-2550CR and SNU-2563CR) by long-term exposure to increasing concentrations of crizotinib, and performed RNA-seq. Both sublines exhibited distinct insensitivities to crizotinib, and a list of differentially expressed genes and their log fold changes together with signaling pathways topology were identified in order to detect the pathways most relevant to the crizotinib resistance. SPIA two-way analysis of RNA-seq revealed multiple signaling aberrations including p53 signaling and focal adhesion (Fig. 4A,B). Specifically, the mRNA level of DKK1 which has been reported to promote the migration and invasion of NSCLC by inhibiting the phosphorylation of β-catenin 19,20 was significantly upregulated in both crizotinib resistant cell lines. The elevated expression of DKK1 was confirmed in mRNA ( Fig. 4C) and protein ( Fig. 4D) level. We further investigated the potential role of increased DKK1 expression in the context of β-catenin signaling. Immunocytochemistry indicated that the nuclear localization of β-catenin was significantly increased in both crizotinibresistant cell lines (Fig. 4E,F).

MPE-derived NSCLC cell lines reveal heterogeneous drug responses caused by molecular diversity.
We assembled an 18-compound library for drug screening. In total, 30 cell lines including two crizotinib resistant sublines from 28 patients were successfully screened in experimental triplicate, generating > 500 measurements of cell line-drug interactions. As a first validation, k-means clustering of AUC values indicated that both cell lines and drugs were likely grouped with their molecular subtypes. For instance, MPE-derived NSCLC cell lines with EGFR mutations such  Table 5). Among the MPE-derived NSCLC cell lines harboring fusion genes, two cell lines with the CD74-ROS fusion gene were selectively sensitive to MK-5108, an Aurora A inhibitor (marked with triangle Fig. 5A). Two crizotinibresistant sublines displayed insensitivities to other next generation ALK inhibitors including alectinib and ceritinib as well. In accordant with the previous study 21 , ceritinib displayed better responses compared to two other ALK-inhibiting reagents. ICG-001, a β-catenin/CBP interaction inhibitor specifically exhibited good response on both crizotinib-resistant sublines (marked with rectangle Fig.s), which recapitulated the role of the increased nuclear localization of β-catenin caused by the elevated DKK1 on the crizotinib resistance.
We further integrated the genetic factors of the MPE-derived NSCLC cell lines with 18 drugs using Wilcox ranked sum test (Fig. 5B). Gene-drug interaction analysis not only reconfirmed that EGFR targeted therapy effectively worked on the EGFR-mutated cell lines, but also revealed that cell lines harboring the EML4-ALK fusion displayed significant insensitivities to EGFR targeted therapy. As indicated earlier, cell lines with CD74-ROS fusion were specifically sensitive to MK-5108 as well as cyclopamine.
Although gene-drug interaction analysis demonstrated the heterogeneous drug responses in accordance with the presence of the molecular target, more comprehensive interpretation was achieved by transcriptomic clustering. We performed multivariate analysis of variance (MANOVA) to find which drugs accorded with the three different transcriptional subtypes. MANOVA indicated that EGFR TKIs were significantly associated with the mRNA expression patterns (Supplementary Table 6). For instance, SNU-2585 and SNU-2588 cell lines were free of EGFR mutations, yet still positioned adjacent to cell lines with EGFR mutations on the transcriptomic analysis. This mRNA expression pattern is reflected by the notably good responses of these cell lines to EGFR TKIs, which suggests that the shape of transcriptomic landscape can be a determinant factor for EGFR-targeting therapy (Fig. 5C). In order to further associate the transcriptomic expression patterns to the drug responses, we then calculated the Pearson correlation coefficient with p values between the transcriptional heterogeneity score (DEPTH) and AUCs of each drug. The AUC value of cyclopamine was in direct proportion to the DEPTH score (Fig. 5D). Then, we integrated the genetic factors with the response to cyclopamine in order to further explain outliers. The mutational status of ERBB2 was highly associated with the correlation coefficient of DEPTH score www.nature.com/scientificreports/ with the AUC of cyclopamine (Fig. 5E). Most of the cell lines with mutated ERBB2 located outside of the confidence interval and had reverse correlation of DEPTH score with the AUC of cyclopamine.

Discussion
Nearly 20% of lung cancer patients develop malignant pleural effusions (MPEs) 22,23 . MPE is an indicator of the poorer survival rate in lung cancer patients (about 7.5 months) compared to those without effusions (about 12-18 months) 23 . In the tumor-node-metastasis (TNM) staging classification of lung cancer, patients with MPE are categorized as category M1a (stage IV) as it is a sign of metastatic dissemination 24 . Although few studies have focused on malignant tumor cells from pleural fluids, only genomic analysis including mutational aspects and copy number variations was conducted [25][26][27] . Giuseppe Roscilli et al. integrated the mutational characteristics of human lung adenocarcinoma cells derived from MPEs to the patients chemosensitivity 6 , yet the number of included drugs and patients were not enough to draw solid conclusion.
In this study, we present 28 MPE-derived lung cancer cell lines and two crizotinib resistant sublines that were comprehensively characterized in accordance with multiple drug responses. These in vitro models retained the representative lung cancer mutations in parallel with mutational profiling of other external datasets consisting of lung cancer cell lines (CCLE) as well as tumor samples (TCGA). Our data revealed that MPE-derived lung cancer cell lines involves higher rates of EGFR mutations and fusion events compared to the lung adenocarcinomaoriginated cell lines. Supporting the notion of distinguished EGFR and fusion aberrations, our transcriptomic analysis demonstrated that MPE-derived cell lines were subtyped in parallel with the genetic aberrations. The transcriptomic subtyping of EGFR-mutant NSCLC has been performed in other studies as well 28,29 . Our study has demonstrated that MPE-derived NSCLC cell lines also retain the transcriptional characteristics of EGFR-mutant NSCLC, which were distinctive to NSCLC cell lines with other driver mutations such as KRAS or fusion genes.
We suggested that the nuclear localization of β-catenin in MPE-derived NSCLC cell lines through increased DKK1 expression potentially caused crizotinib resistance. The acquire crizotinib resistance has been also associated with both MET amplification and exon 14 skipping mutation in lung cancer 30 . We further inspected if MET aberrations occurred in our two crizotinib sublines as well. We confirmed that MET was clearly amplified in SNU-2550CR cell line compared to its parental SNU-2550 cell line regardless of crizotinib treatment, yet no MET amplification was found in SNU-2563 and SNU-2563CR cell lines. This suggested that MET amplification indeed accounted for acquired crizotinib resistance in SNU-2550 set, and the SNU-2563CR cell lines had different resistance mechanism such as nuclear localization of β-catenin. No MET mutation was found in both resistant sublines.
Our work also enabled the integration of highly prevalent genetic and transcriptional characteristics of MPEderived lung cancer cell lines to various drug responses. We re-confirmed that targeted therapies, especially EGFR TKIs showed high sensitivities in accordance with the presence of key mutations using Wilcox ranked sum test (Fig. 5B) and MANOVA (Supplementary Table 6) as reported by several previous studies. Remarkably, www.nature.com/scientificreports/ we identified two EGFR-wild type cell lines (SNU-2585 and SNU-2588) were exceedingly sensitive to EGFRtargeting drugs. Those two cell lines showed even better response to Dacomitinib and Afatinib than EGFR-mutant cell lines (Fig. 5C). Hierarchical clustering analysis using transcriptomic landscape demonstrated that both cell lines were tightly clustered with the EGFR-mutant group (Fig. 3B), which enabled us to conclude that the translational profile can also be a determinant factor for accessing the drug response of targeted therapy. In this perspective, we suggest the downside of using mutation and transcriptomic data separately in clinical diagnosis when predicting MPE-complicated NSCLS patient responses. In analogous with this result, the correlation between the response of MPE-derived NSCLC cell lines to Cyclopamine and transcriptional heterogeneity score (DEPTH) depended on the mutational status of ERBB2 (Fig. 5E).
In conclusion, our work emphasized the significance of integrating the mutational status with transcriptomic subtypes using statistical models such as Pearson correlation coefficient and MANOVA in order for the accurate assessment of drug responses in lung cancer patients with MPE.

Establishment and maintenance of human NSCLC cell lines. The research protocol was reviewed
and approved by the institutional review board of the Seoul National University Hospital (IRB No. 1102-098-357). The study was performed in accordance with the Declaration of Helsinki. Written informed consent was obtained from all patients enrolled in this study.
Cell lines were established from malignant pleural effusion (MPE) derived from pathologically proven nonsmall cell lung cancer (NSCLC). A total of 28 MPE specimens of human NSCLC from 28 different patients who underwent a pleural intervention for the palliation of dyspnea were obtained from Seoul National University Hospital (Seoul, Korea). MPE samples were directly transferred from the operating room to the laboratory for cell culture. Tumor cells were spun down by centrifugation the MPE sample at 300 rpm for 5 min, and re-suspended with Opti-MEM I (GIBCO, CA, USA) supplemented with 5% fetal bovine serum (GE Healthcare Life Sciences, Buckinghamshire, UK) and 1% Penicillin/Streptomycin. Gathered cells were then seeded into T-25cm2 flasks (Corning, NY, USA). Confined-area trypsinization or scraping method was used to attain a pure tumor cell population when stromal cells like mesothelial cells or fibroblasts grew in the initial culture. Established cell lines were sustained in RPMI 1640 medium with 10% fetal bovine serum and 1% (v/v) penicillin and streptomycin (10,000U/ml). Cultures were maintained in humidified incubators at 37 °C in an atmosphere of 5% CO2 and 95% air. The initial passage was assigned when substantial tumor cell growth was detected, and successive passages were given at sub-confluence after trypsinization. When one culture population contains both floating and adherent cells, floating cells were gathered by centrifuging the medium and dispersed by pipetting. All cell lines introduced in this study will be deposited in the Korean Cell Line Bank (http:// cellb ank. snu. ac. kr) at the initial passage in order to be distributed to researchers worldwide.   www.nature.com/scientificreports/ amplified using an AmpFlSTR identifiler Polymerase Chain Reaction (PCR) Amplification Kit (Applied Biosystems, CA, USA). A single cycle of PCR amplified 15 short tandem repeat markers (CSF1PO, D2S1338, D3S1358,  D5S818, D7S820, D8S1179, D13S317, D16S539, D18S51, D19S433, D21S11, FGA, TH01, TPOX and vWA) and an Amelogenin gender-determining marker containing highly polymorphic microsatellite markers. Amplified PCR products were analyzed by an ABI 3500XL Genetic analyzer (Applied Biosystems).

RNA extraction and reverse transcriptase (RT)-PCR.
To obtain RNA, TRIzol (ambion by Invitrogen, CA, USA) was added to cell pellet acquired from cultured lung cancer cell line. After cell lysis was happened, chloroform was added. Then, the sample was vortexed and 4 °C, 12,000 g centrifugation. Each tube containing the sample was divided into aqua phase, interphase, and organic phase. In those aqua phase was transferred into new tube. Isopropanol alcohol was mixed with an equal part of the aqua phase. Mixed sample was centrifuged at 4 °C, 12,000 g after incubated at − 20 °C. RNA was acquired from pellet. The pellet was precipitated by ethanol. RNA was melted in 4th Distilled water containing RNase inhibitor (DEPC). Complementary DNA (cDNA) was synthesized from extracted RNA by reverse transcription kit (Qiagen, MD, USA). Generally measured concentration of RNA was adjusted to 1ug diluted with gDNA wipeout buffer. Extracted RNA concentration was measured nano-spectrometer (Thermo Fisher Scientific, MA, USA). Each 1ug RNA was mixed with RT buffer, RT primer mix, and RTase. Those compounds were heated 42 °C for 30-45 min purposed on annealing, extension. Heated compound was additionally heated 95 °C for 3 min to denature synthetic cDNA. The list of primers we have used are as follows.  TP53 exon4  TCC CCT GCC CTC AAC AAG AT  AAC CTC CGT CAT GTG CTG TG  60   TP53 exon7  CTG GTG GTG CCC TAT GAG CC  AGG AGC TGG TGT TGT TGG GC  61   TP53 exon8  CCT CTT TCC TAG CAC TGC CC  CAA ATG CCC CAA TTG CAG  Sanger sequencing. PCR product was precipitated by 5% sodium acetate buffer (Sigma-Aldrich) and 95% ethanol mixed solution. Then washed product was set on ice for 10 min and centrifuged at 4 °C, 14,000 rpm. Supernatant was discarded and the product was rinsed this time by 70% ethanol and centrifuged 14,000 rpm. Supernatant was discarded then the products were dried using vacuum concentrator (Eppendorf). 10 μL of distilled water was added to dilute precipitated sample. When the product is all diluted in distilled water, cyclic PCR was carried out. Two separate mixtures for forward and reverse sequences were made where they each include 5X sequencing buffer (Applied Biosystems), Big Dye (Applied Biosystems), forward or reverse primer, distilled water, and product from the previous step. Cyclic PCR was carried out with denaturation step at 96 °C, annealing temperature at 55 °C, and elongation at 60 °C for 25 cycles. The cyclic PCR product was then precipitated with www.nature.com/scientificreports/ unprocessed blotting images are provided with supplementary data file. The blotting image for DKK1 was closely cropped near the membrane edges, and we provide replicates performed as well. The MET protein level was confirmed with ChemiDoc Touch Imaging System (Bio-Rad Laboratiores, CA, USA) using consecutive exposing mode with 1-min interval for 5 min. Germany) was used to examine cells. Diverse magnifications were used for different growth patterns and sizes of cells. The intensity of each channel was fixed for the comparison of target protein expression between samples. Digital resolution, scan speed and the number of pictures averaged were set to 1024 × 1024, 40 s per one channel, and 8 pictures respectively. The pictures were focused on the very bottom of the fixed cells for investigating protruding region of cell colonies and the location of β-catenin.
RNA sequencing. Total RNA was isolated from cell lysate using Trizol (Qiagen) and Qiagen RNeasy Kit (Qiagen). Sequencing libraries were prepared using the Illumina TruSeq Stranded Total RNA Library Prep Kit. Fifty-one million reads were obtained from the cell lysates. Following base-calling and alignment with the Tuxedo Suite, rejected reads were analyzed using FusionMap, ChimeraScan, and Defuse with default parameters for RNA and alignment to GRCh37.72. The output was filtered to include in-frame fusions, with at least one rescued read and two unique seed reads, and exclude known, recurrent artifacts.
Whole exome sequencing. SureSelect sequencing libraries were prepared according to manufacturer's instructions (Agilent sureselect all Exon kit 50 Mb) using the Bravo automated liquid handler. Three micrograms of genomic DNA were fragmented to a median size of 150 bp using the Covaris-S2 instrument (Covaris, Woburn, MA). The adapter ligated DNA was amplified by PCR, and the PCR product quality was assessed by capillary electrophoresis (Bioanalyzer, Agilent). The hybridization buffer and DNA blocker mix were incubated for 5 min at 95 °C and then for10 minutes at 65 °C in a thermal cycler. The hybridization mixture was added to the bead suspension and incubated for 30 min at RT while mixing. The beads were washed, and DNA was eluted with 50 ml SureSelect elution buffer (Agilent). The flow cell loaded on HISEQ 2500 sequencing system (Illumina).

Statistical analysis. Statistical analysis was performed using R program version 3.3.1 (R Foundation for
Statistical Computing, Vienna, Austria) with various packages including maftools, PerformanceAnalytics, survminer, survival, iplot, gplot, and lattice. Fisher's exact test was used to analyze GO analysis of various genes. A multivariate analysis of variance (MANOVA) model was applied to the drug response data matrix with various factors such as the mutational status and the three different transcriptional subtypes. Approximate F value, p-value and Pillai's trace score were obtained for each of the factors/drug pairs. A value of p < 0.05 was considered statistically significant. For hierarchical cluster analysis on a set of dissimilarities, each object was assigned to its own cluster, which an algorithm proceeds through iteratively. Two of the most similar clusters are joined at each stage until there is a single cluster. Distances between clusters are recomputed at each stage by the Lance-Williams dissimilarity update formula according to the particular clustering method being used. Clustering methods include: Ward's minimum variance method, complete linkage method, k-means method, and single linkage method.